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The  geomagnetic  cavity  surrounding  the  earth,  called  the  magneto¬ 
sphere,  contains  an  enormous  amount  of  energy  in  the  form  of  electro¬ 
magnetic  waves,  trapped  energetic  particles,  and  the  static  magnetic 
field.  The  waves  and  particles  are  in  a  state  of  dynamic  equilibrium 
as  they  interact  with  each  other  through  various  plasma  instabilities. 

Our  understanding  of  this  equilibrium  state  has  progressed  to  the  point 
where  it  may  be  possible  to  stimulate  wave  amplification  at  ULF  (^0.1-10  Hz) 
and  VLF  (^1-100  kHz)  by  disturbing  the  equilibrium  through  catalytic 
injections  of  plasma  cloudr  or  particle  beams.  Such  a  prospect  has  great 
significance.  First,  it  orfers  a  direct  test  of  the  fundamental  physical 
processes  that  are  believed  to  control  the  state  of  the  magnetosphere. 
Second,  controlled  injections  might  be  used  to  stimulate  artificial 
amplification  of  VLF  or  ULF  emissions  on  command  and  allow  them  co  be  used 
as  a  wide-area  communication  system. 

One  of  the  dominant  plasma  instabilities  in  the  magnetosphere  is  the 
cyclotron-resonance  interaction  tetween  energetic  particles  and  electro¬ 
magnetic  waves.  The  condition  for  the  interaction  to  occur  is  that  the 
motion  of  the  particle  Doppler  shifts  the  wave  frequency  to  its  local 
cyclotron  frequency.  The  net  exchange  of  energy  between  a  band  of  waves 
and  a  group  of  particles  depends  upon  the  shape  of  the  phase-space  dis¬ 
tribution  of  the  particles  and  the  local  frequency  parameters  of  the  plasma 
medium.  The  interaction  is  undoubtedly  responsible  for  significant  wave 
amplification  and  particle  precipitatirn  in  certain  circumstances  and 
appears  to  be  the  process  responsible  for  maintenance  of  an  equilibrium 
configuration  of  waves  and  particles  in  the  magnetosphere. 
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These  qualities  make  this  interaction  a  logical  first  choice  for 
investigation  as  a  means  of  stimulating  artificial  amplification  of  VLF 
and  ULF  waves  in  the  magnetosphere.  The  interaction  is  sensitive  to 
both  the  background  magnetosphere  parameters  and  the  hot  energetic  particle 
distribution.  Modification  of  the  interaction  can  be  achieved  by  the 
introduction  of  localized  clouds  of  plasma  that  reduce  the  phase  velocity 
of  the  waves  and  cause  enhanced  amplification  by  lowering  the  resonance 
velocity.  Alternatively,  direct  injection  of  hot  energetic  particle  beams 
can  also  appreciably  alter  the  local  amplification  characteristics. 

Previous  studies  of  stimulated  amplification  have  been  based  on 
the  linear  perturbation  theory  for  small  wave  amplitudes.  Their  results 
have  indicated  promising  methods  for  generation  of  signals,  but  the 
perturbation  theory  breaks  down  for  realistic  amplification  rates.  Thus 
the  nonlinear  solution  of  the  problem  is  essential  for  a  quantitative 
description  of  useful  power  levels  for  stimulated  emissions.  However, 
the  nonlinear  aspects  of  the  interaction  have  been  particularly  elusive 
because  its  exact  solution  is  described  by  elliptic  functions  in  the 
complex  plane.  This  fact  has  prevented  useful  application  of  the  results 
because  many  side  constraints  and  most  constants  of  the  motion  remained 
undiscovered. 

This  year,  really  significant  theoretical  progress  has  been  achieved 
at  BNW  on  the  exact  nonlinear  problem.  Foremost  was  the  identification  of 
five  new  constants  of  the  motion.  These  constants,  combined  with  two 
previously  known  constants,  complete  the  set  that  is  required  to  fully 
specify  the  behavior  of  the  interaction.  An  entire  new  vista  of  applica¬ 
tions  is  now  accessible.  For  example,  instead  of  integrating  the  plasma 
kinetic  equations  point  by  point  in  the  six-dimensional  phase  space,  it 
is  possible  to  write  formal  solutions  for  the  particle  distribution  that 
convert  the  equations  into  a  straightforward  eigenvalue  problem.  (Imagine 
integrating  the  Schrodinger  wave  equation  over  all  momentum  and  configura¬ 
tion  space  to  ascertain  the  energy  levels  of  hydrogen,  as  opposed  to 
solving  the  much  simpler  eigenvalue  problem.)  Of  course,  the  solution 
of  the  eigenvalue  problem  in  the  complex  plane  is  by  no  means  trivial, 
particularly  for  realistic  conditions  in  the  magnetosphere,  but  a  vital 
step  has  been  accomplished  in  understanding  the  physical  process  and 
evaluating  the  nonlinear  growth  rate. 
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The  mathematical  technique  fc*'  the  nonlinear  solution  is  fully  estab¬ 
lished  in  the  accompanying  document.  A  computer  program  was  also  pre¬ 
pared  to  evaluate  the  closed  form  solution.  The  program  includes  subroutines 
that  probe  all  aspects  of  the  interaction:  particle  trajectories,  arbitrary 
phase  space  distributions,  nonlinear  wave  amplitudes,  and  systematic 
convergence  in  the  eigenvalue  hyperspace.  The  program  contains  several 
alternate  computing,  printing,  and  plotting  subroutines.  Unfortunately 
support  for  the  research  was  exhausted  before  the  program  could  be  fully 
debugged  and  no  numerical  output  is  available  for  magnetospheric  applica¬ 
tions.  (A  few  numerical  solutions  have  been  obtained  for  waves  in  fusion 
plasmas. ) 

The  theory  is  limited  at  present  to  propagation  in  uniform  homogeneous 
media.  For  magnetospheric  applications  the  medium  is  slowly  vaiying  so 
that  numerical  results  from  this  new  nonlinear  solution  would  b*  helpful 
in  the  study  of  stimulated  amplification.  More  significantly,  the  present 
computer  program  is  the  essential  core  for  an  outer  loop  that  links  a 
series  of  cells  which  approximate  nonuniform  and  inhomogeneous  media. 

In  conclusion  our  nonlinear  solution  provides  the  proper  machinery 
for  a  quantitative  assessment  of  various  techniques  for  stimulation  of 
ULF  and  VLF  emissions  in  the  magnetosphere.  We  believe  our  study  should 
be  continued  to  bridge  the  gap  between  the  esoteric  mathematical  expressions 
which  are  now  available  and  quantitative  evaluation  of  specific  plasma 
injection  experiments  in  the  magnetosphere. 


The  research  work  on  this  ONR  Contract  to  Battel le-Northwest  was  performed 
in  part  by  BNW  staff  and  in  part  by  Aerospace  Corp  staff  through  a  subcontract. 
The  information  in  this  report  describes  only  the  BNW  research  which  involved 
$45,000  of  the  program.  The  remaining  $40,000,  less  fee,  supported  the 
Aerospace  research  which  is  reported  in  the  document,  "Wave-Particle  Interactions 
in  the  Magnetospheric  Plasma."  CATR-74 ( 7420) -1 ]  compiled  by  M.  Schulz. 
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ABSTRACT 

The  interaction  between  a  rotary  electromagnetic  plane  wave  propagating 
along  an  axial  magnetic  field  and  an  arbitrary  phase-space  distribution  of 
charged  particles  is  solved  exactly  for  selected  cases.  The  interaction  is 
described  by  the  nonlinearly  coupled  Vlasov-Maxwell  ion-plasma  field  equa¬ 
tions  which  are  integrated  with  the  assistance  of  four  new  constants  of  the 
motion.  The  ion  orbits  in  momentum  space  are  bipolar  doubly  periodic 
eigenfunctions  of  ion  proper  time,  which  are  obtained  in  closed  form  as  the 
difference  between  two  doubly  quasi-periodic  Weierstrass  zeta  functions. 

The  ion  orbits  in  position  space  are  helical-spiral  doubly  quasi-periodic 
functions  of  ion  proper  time,  expressible  simply  in  terms  of  doubly  quasi- 
periodic  Weierstrass  sigma  functions.  The  complete  ion  distributions  are 
flexible  functions  of  six  constants  of  the  ion  motion:  wave-frame  ion 
energy,  transverse  gyro  center,  an  inner  Hamiltonian  correlating  wave-frame 
ion  momentum  with  wave-frame  axial  position,  and  both  first  and  second  axial 
integration  constants.  The  eigenvalue  determination  intricately  interrelates 
the  wave  propagation  vector,  the  wave  amplitude,  the  axial  magnetic  field, 
the  double  periods,  and  the  bipole  separation.  The  solutions  have  potential 
application  to  stimulated  amplification  of  electromagnetic  noise  in  space, 
induced  particle  precipitation  at  the  base  of  the  magnetosphere,  and  radia¬ 
tive  absorption  and  emission  processes  in  fusion  plasmas. 
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INTRODUCTION 


The  exchange  of  energy  between  electromagnetic  waves  and  magnetically 

confined  particles  via  the  cyclotron-resonance  interaction  has  been  the 

object  of  intensive  study  for  many  years1’4  due  to  its  apol icabil itv  in 

space  physics  and  plasma  fusion  research.  The  wave-oarticle 

interaction  is  described  by  the  coupled  Vlasov-Maxwel 1  equations  which 

have  * ell -known  perturbation  solutions  for  small  wave  amplitudes  and 

1  3 

energy  exchange  rates.  However,  in  many  applications,  such  as 

stimulated  wave  amplification  and  particle  precipitation  in  the  macmeto- 

14—18  1213 

sphere  and  radiative  losses  from  fusion  machines,  the  perturbation 

assumptions  are  no  lonqer  valid,  and  the  exact  nonlinear  solution  is 

mandatory. 

Charged  particle  orbits  under  the  influence  of  a  uniform  magnetic 

field  and  circularly  polarized  electromagnetic  wave  are  known  to  be 

1 9 

described  by  complex  elliptic  Integrals.  In  the  wave  frame,  the  particle 

kinetic  energy  and  an  inner  Hamiltonian  were  well  established  as  constants 
20 

of  the  motion,  but  four  other  constants  remained  unknown.  Nevertheless, 

partial  integration  of  the  exact  equations  for  the  wave-particle  interaction 

21  22 

w as  achieved  using  the  available  orbit  properties.  ’ 

In  our  analysis,  all  six  constants  of  the  particle  motion  are  derived 
and  physically  identified.  With  these  relations  it  is  possible  to  formulate 
a  flexible  phase-space  distribution  function  which  satisfies  the  Vlasov 
equation.  Furthermore,  these  constants  also  permit  explicit  integration  of 
Maxwell's  vector  wave  equation  which  reduces  to  a  nontrivial  eigenvalue 
problem  involving  a  multidimensional  complex  hyperspace.  The  method  for 
numerical  evaluation  of  the  proper  set  of  eigenvalues  is  described  and 
exact  solutions  are  given  for  two  and  three  component  plasmas. 
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PLASMA  PHYSICS 


In  the  four-dimensional  space-time  of  special  relativity  spanned  by 
the  line  element 

di^  =  dx^  +  dy^  +  dz^  +  die  t^  =  -c^dx^, 

the  nonlinearly  coupled  Vlasov  and  Maxwell  ion-plasma  field  equations, 

relating  the  electromagnetic  vector  potential  A(fO  and  the  probability 

density  n  (R,P)  for  ions  of  charge  q  and  rest  mass  m  at  four-position 

R  =  (x  y  z  ict)  with  four-momentum  =  m  dR/dx  =  (p  p  p  iE/c),  take 

s  X  y  z 


the  form 


4,13,23-27 


(1) 
(2) 

where  c  is  the  velocity  of  light,  x  is  ion  proper  time,  E  is  ion  energy, 
V*A  is  the  tensor  curl  of  A,  y  is  the  magnetic  permeability  of  free  space 
(in  rationalized  mks  units),  and  the  last  term  sums  the  four-vector 
current-and-charge  sources  over  all  ion  species  s  and  their  four-momentum 
volume  elements  dV(£). 

The  least  intricate  known  method  of  generating  exact  solutions  to 
these  coupled  equations  proceeds  through  two  formidable  steps: 

Step  I.  Postulation  of  an  analytically  tractable  functional  form 
for  the  electromagnetic  vector  potential  A(R),  followed  by  explicit 
integration  of  the  corresponding  orbit  equations  for  each  ion  species, 

dR(t)/dt  =  P(x)/m$,  (3) 

d£(x)/dx  =  qs[VxA(RL)]*P.(x)/ms.  (4) 


[9ns(R,P)/8R]-P/ms  +  [9ns(R,P)/9P]*[qsVxA(R)] .P/ms  =  0, 
-  V.[VxA(R)]  =  pIs/(qsP/ms)ns(R,P)dV(P), 


j 
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Step  II.  Assembly  of  the  corresponding  set  of  solutions  of  the 
Vlasov  equation  (1)  as  a  flexible  function  of  the  constants  of  the  motion 
for  tne  ion  orbits,  followed  by  parameter  adjustment  both  to  satisfy 
Maxwell's  equations  (2)  and  to  meet  positive-definite  probability-density 
requirements. 
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ION  ORBIT  EQUATIONS 

For  a  rotary  electromagnetic  plane  wave  of  vector  potential  amplitude 
a,  propagating  along  a  uniform  axial  magnetic  field  of  magnitude  b,  the 
wave-frame  vector  potential  A(R)  has  the  form 

A(R)  =  [a*cos(Kz)-b*y/2,  a*sin(Kz)+b*x/2,  0,  0].  (5) 

Tne  complex  Lorentz  rotation  relating  the  wave  frame  (unprined  coordinates) 
to  tne  laboratory  frame  (primed  coordinates)  is  given  by 
x  (!  0  0  0  fx] 

y  =  !  I  't,v  •  y\  ,  with  K2  =  (K')2-(W'/c)2,  (6) 


X 

fl 

0 

0 

y 

0 

1 

0 

z 

'  0 

0 

K'/K 

-iw'/cK 

ict 

0 

0 

z  "00  K'/K  iw'/cK  Z*  *  •  V-/W.  w, 

ictj  (0  0  -iw'/cK  K'/K  lict'J 

where  K'  and  w'  are  the  laboratory-frame  propagation  vector  and  angular 

frequency  respectively.  Permitting  the  wave-frame  propagation  vector  K 

to  span  the  full  complex  field  implies  that  the  wave-frame  (z,t) 

coordinates  both  may  be  complex,  although  the  transverse  (x,y) 

coordinates  remain  real. 

A 

Assembling  the  electromagnetic  force-field  tensor  [B,E/ic]  by  taking 

the  tensor  curl  of  the  vector  potential  (5)  yields 

0  b  Ka  sinKz  0’ 

o*A(Rl  =  "b  0  'Ka  cosKz  0  17) 

_Ka  sinKz  Ka  cosKz  o  0  *  [/J 

0  0  0  0( 

A 

The  space-space  components  nere  are  the  magnetic  force-field  tensor  B. 

Tne  space-time  components  are  the  null  electric  force-field  vector  E/1c, 
wnose  vanishing  here  motivates  the  wave-frame  construction.  The  ion  orbit 
equations  (3-4)  thus  become 
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dx/dx  =  px/m, 

dPx/dT  =  (q/m)(Pyb  +  pzKa  sinKz), 

(8) 

dy/dT  =  p  /m, 

dPy/dT  =  -  (q/m) (pxb  +  pzKa  cos  Kz), 

(9) 

dz/dx  =  pz/m, 

dpz/dx  =  (q/m)Ka(-  px  sinKz  +  py  cosKz), 

(10) 

c'ict/dx  =  iE/cm, 

diE/cdx  =  0. 

(ID 

The  ion  species  index  s  may  be  suppressed  (as  here),  without  loss  of 
clarity,  in  any  equation  not  involving  summation  over  ion  species.  It  is 
to  be  noted  that  the  equations  on  the  left  indicate  that  real  transverse 
(x,y)  coordinates  imply  real  transverse  momenta  (p  ,p  )  and  that  complex 
longitudinal  (z,t)  coordinates  imply  complex  axial  momentum  pz  and  complex 
energy  E. 

Constants  of  the  Motion.  The  last  equation  states  that  the  wave- 
frame  ion  energy  E  is  a  constant  of  the  motion.  Multiplying  the 
d(Px»Py»Pz)/dT  equations  by  (Px»Py»Pz)  respectively  and  adding  them  gives 
f/fy.  *  P/Py  +  PzdPz  ■  0.  (12) 

so  that  wave-frame  momentum  magnitude  p  *  /(p^  +  p^  +  p^)  provides  a 

x  y  z 

second  constant  of  the  motion. 

The  structure  of  the  transverse  equations  (8-9)  suggests  merging 
them  in  a  transverse  complex  field,  with  x  +  jy  =  re^  and  p  +  jp  =  pxeJ^, 

y 

to  obtain  as  their  transverse  polar  form, 
d(rej*)/di  =  P^/m,  d(pAeJ>)/dx  =  -  j  (q/m)[bpAe^  +pzKaejKz],  (13) 
where  j  =  /(-l )  has  been  used  for  the  transverse  complex  field  since 
i  =  /(-l)  lies  in  the  longitudinal  complex  field.  The  last  equation  now 
integrates  to 

d[pAe^  +  q(aejKz  +  jbre^)]  =  0,  (14) 

which  shows  that  the  transverse  components  of  the  transverse  quasi-momentum 
pAeJl^  +  q(aejKz  +  jbre^)  provide  a  third  and  fourth  constant  of  the  motion. 


Inner  Hamiltonian.  The  eight-dimensional  orbit  equations  (8-11) 
possess  seven  constants  of  the  motion,  of  which  we  now  have  four.  A 
fifth  constant  of  the  motion,  having  the  character  of  an  inner  Hamiltonian 
coupling  the  intrinsic  physics  of  the  problem,  next  is  obtainable  by 
resolving  the  transverse  momentum  equation  (13)  parallel  and  normal  to 
eJ^  to  give 

dPi/dx  =  (q/m)pzKa  sin(Kz  -  \f>) ,  (15) 

di^/dT  =  - (q/m) [b  +  (p2/p1)Ka  cos(Kz  -  <,)].  (16) 

Expressing  pA  in  terms  of  pz  by  use  of  p2  +  pz  =  p2  converts  the  first 
transverse  equation  (15)  to 

dp2/dr  =  -(q/m)pJLKa  sin(Kz  -  ^).  (17) 

The  coupling  of  the  helical  phase  (Kz  -  )  here  suggests  subtracting 

d(Kz)/dx  from  both  sides  of  equation  (16)  to  obtain 

d(Kz  -t|<)/dx  =  (K/m)pz  +  (q/m)[b  +  (pz/pA)Ka  cos(Kz  -ij>)].  (18) 

The  tra  .sverse  momentum  equations  (17-18)  consist  now  of  two  equations  in 
two  unknowns,  the  helical  phase  (Kz  -  <j>)  and  the  axial  momentum  pz. 

In  this  two-dimensional  subspace,  we  can  use  the  fact  that  any  two- 
dimensional  differential  form  possesses  an  integrating  factor.  Ratioing 
equations  (17-18)  gives  the  differential  form 
(q/m)pAKa  sin(Kz  -ip)  d(Kz  -t|>) 

+  ( ( K/m ) pz  +  (q/m)[b  +  (pz/pA)Ka  cos(Kz  -^)]}  dpz  =  0.  (19) 

Differentiating  the  first  differential  coefficient  with  respect  to  pz, 
using  pA  +  p2  =  p2,  yields 

-(q/m)(pz/pA)Ka  sin(Kz  -^), 

and  mentally  differentiating  the  second  differential  coefficient  with 
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respect  to  (Kz  -  \Jj)  yields  precisely  this  same  result.  The  integrating 
factor  thus  is  unity,  and  the  differential  form  (19)  integrates  directly 
to  provide  our  fifth  constant  of  the  notion,^ 

d[(pz  +  qb/K)2/2m  -  (qa/m)px  cos(Kz  -  i^)]  =  0,  (20) 

which  has  been  scaled  to  the  dimensions  of  energy  to  facilitate  immediate 
physical  interpretation.  The  structure  of  this  inner  energy  constant 
shows  that  when  (K,z,pz)  are  all  pure  real  or  pure  imaginary,  each  ion 
stably  rides  the  rotary  traveling  electromagnetic  wave,  energetically 
trapped  in  its  momentum-dependent  helical  potential  well. 

Further  analysis  is  simplified  considerably  by  shifting  the  four- 
momentum  representation  to  polar  coordinates,  as  defined  by 

E/mc2  =  coshA,  p/mc  =  sinhA,  (21) 

Pz/mc  =  sinhA  cosO,  pA/mc  =  sinhA  sine, 

p  /mc  =  sinhA  sine  cosip, 

A 

Py/mc  =  sinhA  sine  sinip . 

Here  A  is  the  Lorentz  angle  through  which  an  ion  has  been  rotated  from 
rest,  e  is  the  polar  pitch  angle  of  the  ion  momentum  with  respect  t>  the 
axial  magnetic  field,  and  ip  is  the  azimuthal  angle  of  the  ion  momentum 
about  the  axial  magnetic  field.  In  the  wave  frame,  both  Lorentz  ang^  A 
and  pitch  angle  0  may  be  complex,  but  azimuth  angle  is  real. 


4- 

The  existence  of  this  constant  of  the  motion  has  been  noted 
previously  by  Laird2^. 


In  any  two-dimensional  subspace,  it  is  well  known  that  ar,y  constant 

of  tne  motion  can  be  used  as  a  Hamiltonian  by  adroit  choice  of  canonical 

conjugates.  The  polar  momentum  representation  suggests  dividing  the  inner 

energy  (20)  oy  p2/m  to  define  the  inner  Hamiltonian  h  as  the  dimensionless 
ratio, 

and  then 

dimensionless  ratio  form,  to  obtain 

d(cose)/d{Kpx/m)  =  -  (qa/p)sin6  sin(Kz  -  \p),  (23) 

=  -  3h(Kz  -  \p,  cose)/3(Kz  -  ^), 

d(Kz  -  f,)/d(ICpT/m)  =  (cose  +  qb/Kp)  +  (qa/p)cote  cos(Kz  -  *),  (24) 

=  3h(Kz  -  \p,  cose)/3(cose). 

It  Becomes  clear  that  n(Kz  -  \p,  cose)  can  be  used  as  an  inner  Hamiltonian 
function,  with  (Kz  -  ip,  cose,  Kpx/m)  respectively  taken  as  dimensionless 
(position,  momentum,  time)  coordinates. 


Tnis  point  in  our  analysis  logic  marks  the  outer  boundary  of 
previously  existing  firm  theory,  as  delineated  by  our  literature  search 
for  recent  related  work  referenced  above  as  well  as  others2**"'^. 


h  =  (cos0  +  qb/Kp)2/2  -  (qa/p)sin0  cos(Kz  -  ^) , 
rescaling  the  transverse  momentum  equations  (17-18)  to 
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FIRST  INTEGRALS 

Axial  Momentum.  The  inner  Hamiltonian  equations  of  motion  (23-24) 
consist  of  two  equations  involving  two  unknowns,  the  helical  phase  (Kz  -4)) 
and  the  wave-frame  pitch  cosine  (cosG).  Maneuvering  toward  one  equation 
in  one  unknown  leads  to  eliminating  tne  helical  phase  (Kz  -  ij;)  by 
combining  the  Hamiltonian  constant  (22)  with  the  square  of  the  pitch 
cosine  derivative  (23),  to  obtain 

[d(cose)/d(KpT/m)]2  =  (qa/p) 2( s i n2e ) [  1  -  cos2(Kz  -  1^)], 

=  (qa/p)2(l  -  cos20)  -  [(cose  +  qb/Kp)2/2  -  h]2. 

We  now  nave  one  equation  involving  one  unknown,  the  pitch  angle  cosine 
(cose).  Simplifying  the  most  complicated  term  with  an  origin  shift  yields, 
sequential ly, 

u  i  cose  +  qb/Kp,  (25) 

[du/d(Kpi/m)]2  =  (qa/p)2[l  -  (u  -  qb/Kp)2]  -  (u2/2  -  h)2, 

[du/d(Kpt/m)]2  =  -  u4/4  +  [h  -  (qa/p)2]u2  +  2(qa/p)2(qb/Kp)u 

+  (qa/p)2[l  -  (qb/Kp)2]  -  h2.  (26) 

The  square  of  the  derivative  of  u  on  the  left  here,  with  the  fourth  degree 

polynomial  in  u  on  the  right,  leads  to  immediate  recognition  that  this 

equation  is  soluole  in  closed  form  in  terms  of  doubly  periodic  Weierstrass 
37-41 

functions. 


* 


Maneuvering  toward  the  Weierstrass  differential  equation, 

[dP(z)/dz]2  =  4(P-ei)(P-e2)(P-e3),  e^e^  =  0,  P(z=0)  =  1/z2,  (27) 

by  first  factoring  the  quartic  polynomial  on  the  right  of  equation  (26) 
and  then  factoring  out  one  of  the  factors,  yields  the  reduction 

[du/d(iKpx/2m)]2  =  (u-u] ) (u-u2)(u-u3) (u-u4) ,  (28) 

[-d(u-u4)_1/d(iKpT/2m)]2  =  [(u-u1 )/(u-u4)][(u-u2)/(u-u4)][(u-u3)/(u-u4)] , 

=  [1-(u1-u4)/(u-u4)][1-(u2-u4)/(u-u4)][1-(u3-u4)/(u-u4}]. 

The  cubic  polynomial  structure  here  involving  (u-u4)_1  suggests  scaling  to 
w  =  -(u1-u4)(u2-u4)(u3-u4)/4(u-u4),  to  obtain  (29) 

[dw/d(iKpT/2m)]2  =  4(w-w-|  )(w-w2)(w-w3) ,  where  (30) 

w1  =  -(u2-u4)(u3-u4)/4,  w2  =  -(u3-u4)(u1-u4)/4,  w3  =  -(u^) (u2~u4)/4.  (31) 

An  origin  shift  to  the  centroid  of  the  cubic  roots, 

w  =  (w1+w2+w3)/3,  e1  =  wrw,  e2  =  w2-w,  e3  =  w3~w,  (32) 

now  yields  the  Weierstrass  equation  (27)  with 

P(iKpx/2m+X)  =  w-w,  (33; 

where  X  is  a  flexible  integration  constant.  Backtracking  through  the 
preceding  transformation  hierarchy  to  assemble  the  first  integral  for 
the  ion  orbits  yields,  finally, 

cosG  +  qb/Kp  =  u4  -  (u1-u4)(u2-u4)(u3-u4)/4[P( iKpx/2m+X)  +  w].  (34) 

The  quartic  roots  (u^ ,u2>u3,u4)  obviously  can  be  related  to  the 
coefficients  of  the  quartic  polynomial  on  the  right  of  equation  (26)  by 
the  symmetric  polynomial  forms 
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(35) 

(36) 

(37) 

(38) 

Parameter  Inversion.  The  powerful  general  theory  of  doubly  periodic 
functions  can  be  used  to  clean  up  the  curiously  asymmetric  and 
computationally  awkward  parameterization  appearing  in  the  solution  (34). 

The  structure  of  this  doubly  periodic  so  ution  displays  two  first  order 
poles  in  each  primitive  period  cell,  situated  at  the  two  (and  only  two) 
primitive  roots  of 

P(6)  =  -  w,  which  lie  at  iKpi/2m  +  X  =  ±  6.  (39) 

Combining  (31,32,39)  with  (27)  yields  the  parametric  interrelationship, 
-[(u1-u4)(Uj,-u4)(u3-u4)]Z/43  »  Ce1-P(6)3[e2-P(«)][e3-P(fl)D . 

-  -[P'(6)]2/4,  or 

(u1-u4)(u2-u4)(u3-u4)/4  -  -P 1 (6) .  (40) 

The  choice  of  sign  in  the  last  step  completes  the  definition  of  the  bipolar 
separation  (2fi)  by  sharpening  the  (♦)  choice  left  by  (39).  Our  solution 
(34)  thus  can  be  parameterized  more  neatly  now  as 

cosO  +  qb/Kp  u4  +  P'  ( <5 )/ [P ( iKpi/2m+A)  -  P( 6 ) ] ,  (41) 

with  6  parametrically  replacing  both  w  and  (uj-u4)(u2-u4)(u3~u4)  via 
equations  (39-40). 

Shifting  the  remaining  parameterization  to  the  symmetric  Weierstrass 
parameters  (e^.e^e^),  in  an  effort  to  get  rid  of  the  quartic  root  u4 
surviving  as  an  obnoxious  asymmetry  in  the  solution  (41),  leads  to  combining 


equations  (31,32,39)  to  obtain 


erP(5)  =  -(u2-u4)(u3-u4)/4, 
e2-P(5)  =  -(u3-u4)(u1-u4)/4, 
e3-P(6)  =  -(u1-u4)(u2-u4)/4, 
and  then  making  use  of  equation  (40)  to  extract 


Tu4  =  P'(6)/[erP(6)], 

(42) 

2-u4  =  P'(6)/Ce2-P(6)], 

(43) 

3-u4  =  P'(5)/[e3-P(6)]. 

(44) 

Adding  these  three  equations,  using  uj+u2+u3+u4  =  0,  yields  finally, 

u4  =  {P,(6)/[P(6)-e1]+Pl(6)/[P(6)-e2]+P'(5)/[P(6)-e3]}/4, 

-  (In  [P'(6)]2}'/4,  from  equation  (27);  hence 


u4  =  P" ( 6)/2P 1 (6) , 

(45) 

and  from  (42-44)  then, 

u1  =  FU(6)/2P'(6) 

-  P'(6)/[P(«)-e,], 

(46) 

u2  =  P*'(6)/2P'(6) 

-  P'(S)/[P(S)-e2], 

(47) 

u3  =  P" ( 6 )/2P ' (6) 

-  P'(6)/[P(6)-e3J. 

(48) 

Double  differentiation  of  the  Weierstrass  half-period  addition 
theorem^(0-79)  ,  with  implicit  cyclic  (1  2  3)  index  permutation,  yields, 
sequentially, 


[P(z+ui1)-e1][P(e)-e1]  =  (e^  Me-j-ej ), 

P,(z+w1)/[P(2+w1)-e1]+P'(z)/[P(z)-e1]  =  0, 

P"(a+w^ )/ P 1 (z+wi )  =  P"(z)/P,(z)-2P,(2)/[P(z)-e1].  (49) 

Reference  38  (Equation  0.79). 
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Here  (<V“2.<‘.3)  are  the  Weierstrass  half  periods,  with 

P(“')  =  el-  p(“2)  =  e2,  P(o)3)  =  e3, 

THe  quartic  roots  (-,,-48)  thus  reduce  neatly  to 


u,  = 


u„  = 


u„  = 


Un  = 


P"(6+tJ])/2P1(6+a)1), 
Pn(6+a32)/2P'(6+a>2), 
P"(  6+^3  )/2P '  ( 6+u>3  ) , 
P"(6)/2P' (6) , 


(50) 

(51) 

(52) 

(53) 


and  our  bipolar  solution  (41)  now  becomes 

cosO  +  qb/Kp  =  P"(4)/2p.(6)  +  P-(5)/[P(lKpT/2m+A)  _  p(J)]_ 

Further  integration  is  facilitated  by  shifting  this  representation 
o  He  Weierstrass  additive, y  quasi-periodic  zeta  function,  as  defined  by 
dc(z)/dz  -  -  PCr),  ?(e.0)  a  1/e. 
by  use  of  the  transparent  identity  38(°-9°82) 

»'(*.)/[PW-P(..)].t(i.et).t(Wl|Jt 
which  yields 

COSO  +  qb/Kp  «  4(iKPT/2m+A-o)  .  C< iKPT/2m+A+5)  +  2?(d)  +  p„({)/2p,. 

Use  of  the  further  identity33^.*. 7)  ' "  U) ' 

He)/2P'W  «  C(2e)  -  2c(e), 

reduces  our  bipolar  doubly  periodic  solution  to  its  most  useful  form, 

.  7  '  ParaWteri2ed  f°r  b°‘h  investigation  and  numeric 

evaluation. 


cosO  + 


qWKp  =  «UKpT/2mH-d)  -  c(iKpT/2m+A+6)  7  ?(^7 


(56) 
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This  key  first  integral,  in  brief  resume',  gives  the  wave-frame 
pitch-angle  cosine  (cose)  at  ion  proper  time  t  as  a  bipolar  doubly 
periodic  function  exp -essed  explicitly  as  a  phase-shift  difference  of  the 
doubly  quasi-periodic  Weierstrass  zeta  function.  The  primitive  (tiniest) 
periods  (2^  ,2^)  and  the  primitive  (tiniest)  bipole  separation  (26)  are 
parametrically  related  to  the  three  coefficients  (h,qa/p,qb/Kp)  appearing 
in  the  differential  equation  of  motion  (26)  by  the  cascade  of  boxed 
equations  (35-38,50-53). 

Transverse  Momentum.  Transverse  momentum-projection  first  integrals 
now  can  be  generated  from  the  inner  Hamiltonian  (22)  and  the  transverse 
equation  of  motion  (23),  respectively,  as  the  rotary  orthogonal  projection 
cosines , 

sinG  cos(Kz  -  ip)  =  (p/qa)[(cos0  +  qb/Kp)2/2  -  h],  (57) 

sine  sin(Kz  -  ip)  =  -  (p/qa)  d(cos0)/d(Kpx/m) .  (58) 

Forming  the  square  of  the  axial  first  integral  (56)  for  subsequent 
insertion  in  (5  ),  by  Laurent  expansion  about  either  pole,  yields  as 
leading  terms, 

cose  +  qb/Kp  =*  ±l/(iKpx/2rn+A?6)  ±  P(26) •  (iKpx/2m+A+6) , 

(cose  +  qb/Kp)2  =  l/(iKpx/2m+A+6)2  +  2P(26). 

From  the  general  theory  of  doubly  periodic  functions,  it  follows  that 

(cosO  +  qb/Kp)  =  P(iKpx/2m+A-6)  +  P(iKpx/2m+A+6)  +  P(26).  (59) 

Differentiating  the  axial  first  integral  (56)  for  subsequent  insertion  in 
(58)  yields  at  once, 

d(cos6)/d(Kpi/m)  =  -  ( i/2 ) [P ( i Kpx/2m+A-6 )  -  P( iKpx/2m+A+6)] .  (60) 


Assembling  these  pieces  now  gives  the  transverse  momentum-projection 
first  integrals  (57-58)  as  the  bipolar  doubly  periodic  functions  of  ion 
proper  time  x. 


(61) 
(62) 

The  astonishing  structural  similarity  of  these  rotary  orthogonal 
projectors  suggests  multiplying  the  second  by  i  and  then  adding  and 
subtracting  the  two  to  obtain  the  cisoidal  resolution, 

(63) 

(64) 


sine  e'*Kz"**  *  (p/qa)[P(iKpT/2m+X+6)+P(26)/2-h], 
sine  e‘1(Kz"4')  =  (p/qa)[P(iKpT/2m+A-6)+P(26)/2-h]. 


sinG  cos(Kz-i^)  =  (p/2qa)[P(iKpx/2m+A-6)+P(iKpx/2m+A+6)+P(26)-2h], 
sinG  sin(Kz-i/i)  =  i (p/2qa)[P(iKpx/2m+A-6)-P(iKpx/2m+A+6)] . 


Ratioing  this  cisoidal  pair  extracts  the  helical-phase  first  integral  as 


to  display  the  explicit  solution  for  the  second  inner  Hamiltonian 
equation  of  motion  (24). 
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SECOND  INTEGRALS 

ftxjA1  Position.  The  first  of  the  axial  equations  of  motion  (10) 
combines  with  the  axial  first  integral  (56)  to  give 

d(Kz)/d(KpT/m)+qb/Kp  =  c(iKpx/2m+A-6)-c(iKpx/2m+A+6)+i;(25). 
Integrating  by  use  of  the  mul tipi icatively  quasi-periodic  Weierstrass 
sigma  function,  as  defined  by 

(66) 


(67) 

This  axial  second  integral  gives  the  wave-frame  ion  axial  position  z(t) 
as  an  explicit  doubly  quasi-periodic  function  of  ion  proper  time  t,  wi.th 
z0  =  z(0). 

From  the  multiplicative  double  quasi-periodicity  of  the  Weierstrass 
sigma  function38^0’73) 


[da(z)]/a(?)d?  =  £(?),  a(z^0)  ^  ?,  yields 


e’'K <^“)  ■  e[t(2S)-^^(iKpx/m).o2(1KpT/2m+A.{)o2a+6) 


/cT(iKpx/2m+A+6)a  (A-6) 


oU+v^j+v^)  =  -  e2['>,s(^1  )+v2;(<o2)]a  (68) 

with  sPanning  all  positive  and  negative  integers  except  the  (0,0) 

pair,  it  follows  that  our  axial  position  integral  (67)  has  the 
multiplicative  double  quasi-periodicity 
e  i  K[  z  ( i  Kpi/2m+v^w^  +v2w2  )  -z  ( i  Kpx/2m- -v^ )  ] 

=  e4^^26)'^b/KP^viu,i+v2t02^8^vl^(a)l)+v2;^2^6  (69) 

The  helical  concept  of  axial  pitch  thus  can  be  ascribed  to  this  discrete 


doubly-indexed  subset  of  ion  orbits,  as  given  explicitly  by 
z  ( i  Kpx/2m+v1  w1  +v2w2 )  -z  ( i  Kpx/2m- -v^ ) 

=  ■f4^(26)-9b/|<P](v1a)1+v2u)2)-8[v1;(u)1)+v2;(w2)]6}/iK. 


(70) 
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This  remarkable  equation  states  that,  as  the  ion  momentum  cycles  through 
any  doubly-indexed  Weierstrass  period  2v1w^+2v2(jJ2»  the  wave-frame  ion 
axial  position  z  changes  by  the  fixed  axial  pitch  appearing  here  on  the 
right.  A  fascinating  observation  of  fundamental  relevance  is  that  all 
ion  orbits  accessible  to  numeric  evaluation  on  a  digital  computer  are 
contained  in  this  discrete  doubly-indexed  subset,  since  digital  computer 
numbers  span  only  discrete  integer  ratios. 

Transverse  Position.  The  transverse  second  integrals  now  can  be 
assembled  quite  simply,  by  first  pulling  the  transverse  qua si -momentum 
constants  of  the  motion  (14)  out  of  the  transverse  complex  field  in  the 
(x,y)  component  form  as 

p  sinO  cosij>  +  qa  cosKz  +  qbr  cos(<j>+rr/2), 
p  sin0  si n^  +  qa  sinKz  +  qbr  sin (4>+tt/2) , 
and  then  recombining  them  in  cisoidal  form  in  the  longitudinal  complex 
field  as 

p  sin0  e1'^  +  qa  eiKz  +  qbr  ei  (4>+^/2)^ 
p  sine  e_i^  +  qa  e~lKz  +  qbr 

Keeping  in  mind  here  that  (p,6,K,z)  may  be  complex,  and  making  use  of  the 
axial  first  and  second  integrals  (63,64,67),  yields  for  the  transverse 
position  integrals, 

reH  =  roe1*0  +  i(a/b)(eiKz-eiKz°) 

+  i(p/qb){[sin0  e"1^Kz"'|J^elKz-(p/qa)[P(X-5)+P(25)/2-h]elKZo},  (71) 

re"i<})  =  r0e"i(^0  -  i(a/b)(e"lKz-e"lKz°) 

-  i(p/qb){[sin0  e1  (Kz'^]e"lKz-(p/qa)[P(X+5)+P(26)/2-h]e"lKz°},  (72) 
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where  the  integration  constants  r0,<f>0,z0  refer  to  proper  time  t  =  0  for 
each  ion.  The  factors  rn  the  right  [sine  and  e11^2]  refer  to 

the  explicit  functions  of  ion  proper  time  x  qiven  by  the  axial  first  and 
second  integrals  (63,64,67).  The  periodicity  of  the  former  and  the  aucsi- 
periodicity  of  the  latter  suqgest  gathering  all  the  constant  terms  to  the 
left  side  of  equations  (71-72),  as  a  complex  gyration  center,  and  shifting 
to  complex  polar  coordinates  (R , )  about  this  qyration  center,  as  defined 
by 

Re1*  =  re14>-r0e14>0+i(a/b){l+(p/qa)2[P(x-<5)+P(26)/2-h]}elKz°,  (73) 

Re-1*  =  re"1$-r0e"Ho-i(a/b){l  +  (p/qaf[P(x+6)+P(20/2-h]}e‘lKz°,  (74) 

to  reveal  the  functional  structure 

Re1"  =  i  (a/b ) [1  +  (p/qa)  sine  e”1  ^]elKz, 

Re-1*  =  -i  (a/b)[l  +  (p/qa)  sinO  ej  (Kz-^]e‘iKz. 

Explicit  insertion  of  the  axial  first  integrals  (63,64)  for  sine  e±i(KZ"'J») 
converts  these  transverse  second  integrals  to  the  more  readily  interpre¬ 
table  quasi-helical  form, 

(75) 

(76) 

The  helical  geometric  structure  (R,<i>,z)  on  the  left  combined  with 
the  bipolar  doubly  periodic  proper  time  (x)  structure  on  the  right  show 
that  the  ion  orbits  do  have  the  primary  structure  of  helices  but 
also  have  multi-lobe  oscillatory  secondary  spiral  structure  superimposed 
on  the  heMx.  The  actual  physical  shape  of  these  double  helix  orbits 
depends,  of  course,  on  the  relative  magnitudes  of  the  constants. 


Rei(*-Kz)  =  i(a/b){1+(p/qa)2[p(iKpT/2m+x-5)+P(25)/2-h]}, 
Re-iO-Kz)  =  -i(a/b){l  +  (p/qa)2[P(iKpi/2m+X+6)+P(26)/2-h]} 
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The  striking  similarity  of  these  helical  spiral  ion  orbits  to  the 
spiral  staircase  quantum  structure  of  biophysics  suggests  that  the 
presence  of  a  planetary  magnetic  field  may  play  an  essential  role  in  the 
evolution  and  reproduction  of  life,  by  providing  a  magnetic  axis  about 
which  ions  tend  to  cluster  in  a  spiral  staircase. 

Forming  the  product  and  ratio  of  the  transverse  second  integrals 
(75,76)  extracts  both  the  transverse  gyration  radius  R  and  its  helical 
azimuth  $-Kz  as  the  explicit  multi-polar  doubly  periodic  functions  of 
ion  proper  time  x, 

R2  =  (a/b)2{l+(p/qa)2[P(iKpx/2m+X-6)+P(26)/2-h]} 

*{l  +  {p/qa)2[P(i  Kpx/2m+X+6)+P(26)/2-h]} ,  (77) 

ei2(*-Kz)  a  .{i  +  (p/qa)2[p(iKpT/2m+X-6)+P(26)/2-h]} 

/{l+(p/qa)2[P(iKpx/2m+X+6)+P(26)/2-h]}.  (78) 


ION  DISTRIBUTION 


Viewing  the  eight-dimensional  Vlasov  differential  operator  appearing 
in  equation  (1)  as  simply  the  ion  proper  time  derivative, 

dns(R(x),Pjx))/dx  =  0,  (79) 

makes  it  clear  that  its  complete  solution  consists  of  any  function  of  any 
seven  functionally  independent  constants  of  the  ion-orbit  motion.  Five 
constants  of  the  ion  motion  have  been  displayed  explicitly  in  (11,12,14,22) 
as  [E,p,p,c '^+q(aeJ  +jbre^),h],  and  two  more  have  emerged  implicitly 
in  (56,6/)  as  the  first  and  second  axial  integration  constants  (A,z0). 
Assembling  all  seven  constants  of  the  ioi,  motion  here  for  subsequent 
reference,  as  explicit  or  implicit  functions  of  the  ion  phase-space 
coordinates  (x,y,z,ct;mc,A,0,i|j) ,  yields  in  our  complex  wave  frame: 


(80) 

(81) 


(82) 

(83) 


(84) 


(85) 

(86) 


Constants  of  Ion  Motion  (Wave  Frame) 

i°H  Energy  (E):  From  (11,21),  E  =  me2  coshA. 

jon  Momentum  Magnitude  (p):  From  (12,21),  p  e  me  sinhA. 

ISSl  Ir.a ns verse  Gyro  Center  (XQ,Y0) :  From  (14,21), 

x9  E  x  +  (mc/qb)  sinhA  sin©  sini|>  +  (a/b)  sinKz, 

Yq  h  y  -  (mc/qb)  sinhA  sin9  cos^  -  (a/b)  cosKz. 

Ion  Inner  Hamiltonian  (h):  From  (21,22), 

h  -  [cosO  +  qb/(Kmc  sinhA)]2/2  -  [qa/(mc  sinhA)]  sine  cos(Kz-t^). 

Fjrst  and  Second  Axial  Integration  Constants  (X,z0): 

From  (21,56,67),  using  left  (11)  integral  Kpx/m  =  Kct  tanhA, 

cosO  +  qb/(Kmc  sinhA)  e  c((iKct/2)tanhA+X-6) 

-  S((iKct/2)tanhA+A+6)  +  c(26), 
eiKz0  E  eiKz-[c(26)-qb/(Kmc  sinhA)]iKct  tanhA 

2 

•O  ( (iKct/2)tanhA+A+<$)*o2(A-6)/a2((iKct/2)tanhA+A-6)»o2(A+6) . 


The  first  axial  integral  (85)  gives  the  first  axial  integration  constant 
as  an  implicit  function  X  =  x(ct,mc,A,e).  The  second  axial  inteqral  (86) 
then  gives  the  second  axial  integration  constant  as  a  mixed  explicit-implicit 
function  z0  =  z0(z,ct,mc,A,x(ct,mc,A,e)).  Forming  the  accessible  momentum 
phase  element  from  (21,80,81)  as 

3  2 

dp  dp  dp  dE/c  =  6(m'c-mc)  (m'c)  d(m'c)  sinh  AdA  sinede  dij;,  (87) 
x  y  z 

now  yields  as  the  complete  probability  density  element  for  the  Vlasov  ion 
plasma , 

(88) 

where  the  flexible  ion  nhase-space  distribution  n(X. ,Y. ,z0,A,h,x)  spans  the 

u  o 

six  constants  of  the  motion  remaining  after  meetinq  the  delta-function  rest- 
mass  constraint  in  (87). 

Development  Sequence.  The  establishment  of  a  qeneral  distribution  in 
incorporating  all  of  these  available  constants  is  a  formidable  task.  Certain 
dependences  are  ignored  initially.  For  example,  it  is  important  to  recognize 
that  the  dependence  of  the  ion  distribution  upon  the  gyro-center  (Xn,YQ)  con- 
stants  (82-83)  implies  transverse  (x,y)  variation  of  the  ion  density,  and 
that  dependence  of  the  ion  distribution  upon  the  implicit  first  or  second 
axial  integration  constants  (x,z0)  of  (85-86)  implies  intractably  intricate 
longitudinal  (ct,[z,ct])  variation  of  the  wave-frame  ion  densitv.  The  com¬ 
plementary  inference  here,  then,  is  that  the  transversely  uniform  subgroup 
of  rotational  modes  which  stably  ride  the  rotary  electromagnetic  wave  with 
no  wave-frame  time  variation  should  be  given  by  an  ion  distribution  depen¬ 
dent  upon  only  Lorentz  angle  A  and  inner  Hamiltonian  h,  of  the  sort 

O 

d  N/dxdydz  =  n(A,h)  sinh  AdA  sinede  d^ . 


(89) 


Realistic  further  theoretical  investigation  should  perhaps  focus  first 
upon  this  latter  form,  as  a  fundamental  mode  apt  to  be  of  dominant 
importance  in  most  physical  applications. 

Subsequent  generalization  might  logically  focus  next  upon  transversely 
variant  modes  dependent  also  on  the  gyro-center  (Xfi,Y0)  constants,  of  the 
sort 

d6Ns/dx'f>dz  =  ns(^Q » Yq , A , h )  sinh2AdA  sinede  dip. 

Transversely  variant  modes  like  this  quite  certainly  are  needed  to 
investigate  fusion  containment  leakage  through  the  lateral  wall  of  a 
magnetic  bottle. 

As  a  later  development  phase,  substantial  improvements  in  ion-plasma 
perturoation  tneory  snould  be  obtainable  in  tractable  form  from  the 
complete  distribution  (88),  with  direct  inclusion  of  the  intricately 
coupled  stability  effects  of  perturbing  all  the  constants  of  the  ion 
motion  (X0,Y0,zo.A,h,X). 

Finally,  nothing  more  formidable  than  cascaded  implicit  transformation 
glgeura  seems  to  oe  needed  to  fit  the  complete  distribution  (88)  directly 
to  an  arbitrary  laooratory-frame  initial  distribution,  with  subsequent 
generation  of  tne  full  time-dependent  evolution  of  the  laboratory-frame 
ion  distribution  in  closed  form. 


ELECTROMAGNETIC  FIELD  EQUATIONS 

Ion-Orbit  Current  and  Charge  Integrals.  For  the  axially-guided 
rotary-wave  vector  potential  (5)  and  the  corresponding  complete  Vlasov 
ion  distribution  (88),  Maxwell's  equations  (2)  become 


K£a  cosKz  =  ycj  q  Jn  (Xfi,Yfl,z0 ,A,h,x)  sinh^AdA  sin29de  cosifftty. 


-s'V  V'V'e 


K^a  sinKz  =  PcJsqs/ns(X0,Y0,za ,A,h,X)  sinh^AdA  sin20d0  simj>dij>, 
0  =  VjCIscls/ns(X0»Y0»zo»A.h,X)  sinhJAdA  cocB  sin0de  dip , 

0  =  lJcIsqs/ns(XQ,YQ,Zo ,A,h,A)  coshA  sinfTAdA  sin0d9  dip. 


(91) 

(92) 

(93) 

(94) 


Here  the  ion  distribution  function  ns(X0,YQ,zo  ,A,h,X)  is  complete,  the 
functions  (X@,YQ,h)  are  given  explicitly  by  (82-84),  the  functions  (X,z0) 
are  given  implicitly  by  (85-86),  the  complex  integration  domain  for 
(A,0)  is  obtainable  by  backtracking  to  the  laboratory  frame  via  (6), 
the  azimuthal  (\p)  integration  spans  (-tKi/kti),  and  the  summation  on  s 
spans  all  ion  species  present.  Resolving  the  transverse  (x,y)  equations 
(91-92)  parallel  and  normal  to  the  rotary  vector  potential,  by  multiplying 
the  pair  from  the  left  by  the  rotation  matrix 
f  cosKz  sinKz) 


-sinKz  cosKz 


gives  the  simpler  transverse  pair, 


K  3  *  UcIsqs/ns(X@,Y0  »ze »A,h ,x)  sinh3AdA  sin20d0  cos (Kz-ip)dip,  (95) 


0  =  ">jCXsqs/ns (X0, YQ,ze ,A,h,X)  sinh3AdA  sin20d0  sin(Kz-ip)diJ>.  (96) 


All  four  equations  (93-96)  now  can  be  integrated  explicitly  by 
shifting  to  ion  orbit  coordinates  with  the  differential  transformation, 


sinede  di^  =  d(Kz-tp)  d(cos8)  =  1 3(Kz-<Kcos9)/3(Kpx/m,h) |  d(KpT/m)  dh. 

The  transformation  determinant  here  is  unity,  for  the  inner  Hamiltonian 
equations  (23,24)  evaluate  it  quite  directly  as 

[3  (Kz-tp)/3  (Kpx/m)][3(cos9)/3h]-[3  (cos9)/3(Kpx/m)][3(Kz-i|>)/3h] 

=  [3h/3(cosG)][3(cos9)/3h]+[3h/3(Kz^)][3(Kz-<|>)/3h]  =  3h/3h  =  1. 

The  differential  transform  to  ion  orbit  coordinates  thus  becomes  simply 
sin9d9  dip  =  d(Kz-i^)  d(cos0)  =  d(Kpi/m)  dh,  (97) 

and  making  use  of  the  transverse  and  axial  first  integrals  (61-62,56)  as 
they  arise  now  converts  Maxwell's  equations  (95-96,93-94)  to  their 
elegant  ion-orbit  integral  form, 

K2a  =  (uc2/2a)JsmsJns(X8,Y9,zo,A,h,A)[P(iKpi/2m+X-6) 

+P(iKpx/2m+X+6)+P(26)-2h]  d(Kpi/m)  dh  sinh\dA,  (98) 

0  =  -i(nc2/2a));smsJns(X@,Ye,Zo,A,h,X)[P(iKpT/2m+X-6) 

-P(iKpx/2m+X+6)]  d(Kpx/m)  dh  sinh4AdA,  (99) 

0  “  lJCAscls/ns(Xg,Yg,Zo,A,h,X)[c(26)-q  b/(Km  c  sinhA) 

+C(iKpi/2m+X-6)-^(iKpi/2m+X+6)]  d(Kpx/m)  dh  sinh^AdA,  (100) 

0  ~  ^JCIs(l<;/ns(Xg,Yg,Zo,A,h,X)  d(Kpx/m)  dh  coshA  sinh2AdA.  (101) 

From  the  doubly  periodic  structure  of  both  the  axial  and  transverse 
first  integrals  (56,61-62),  it  follows  that  the  ion  proper  time  integrals 
here  must  span  one  doubly-indexed  Weierstrass  period  2v1w1+2v2io2  ,  to 
properly  span  the  momentum  angular  phase  space  once  and  only  once. 

The  transverse  current  integral  (99)  normal  to  the  rotary  vector 


24 


potential  manifestly  vanishes,  by  integrating  with  the  help  of  the  zeta 
function  (55)  to  a  doubly  periodic  zeta  function  difference  having  the 
same  value  at  both  ends  of  the  Weierstrass  period.  With  the  help  of  the 
zeta  function  (55)  and  its  additive  double  quasi-periodicity32^*^3^, 

C (z+ViU)i  +\>2^2 )  =  C ( 2— v-jto-j 2\>i^(to-j)  +  2\>2 C (^>2 ^ *  (102) 

and  with  the  help  of  the  sigma  function  (66)  and  its  multiplicative  double 
quasi-periodicity  (68),  Maxwell's  remaining  three  equations  (98,100,101) 
integrate  mentally  now  to  the  merged  Vlasov-Maxwell  ion-plasma  triad. 


(103) 

(104) 

(105) 

All  quantities  appearing  in  these  equations,  except  (K,a,b,y,c), 
implicitly  contain  a  suppressed  ion-species  label  s.  The  first  of  this 
triad  (103-105)  relates  the  rotary  electromagnetic  wave  amplitude  to  the 
rotary  ion  current  source.  The  last  two  state,  respectively,  that  both 
the  axial  ion  current  and  the  charge  density  are  zero. 

The  ultra-intricate  dependence  of  the  Weierstrass  primitive  periods 
(2uj-j , 2oj2 )  and  of  the  primitive  bipole  separation  (26)  upon  the  three 
coefficients  [h,qa/(mc  sinhA),qb/(Kmc  sinhA)]  in  the  axial  equation  of 
motion  (26),  via  the  cascade  of  boxed  equations  (35-38,50-53),  makes 
further  integration  at  this  full  level  of  completeness  quite  formidable. 


iK2a2/2y  =  ^smsc2/ns(Xe,Yg,Zo,A,h,A){[P(26)-2h](v^w-|+v2(jJ2) 
-2[v15(w1  )+o2c(a)2)]}  dh  sinh^AdA, 

0  =  ^sM^V  Vz°’A,h’^{^(26)-qsb/(Kmsc  sinhA)](v1<Vv2(JJ2) 
-2[v-|  c(w-|  )+v2?(w2)]5)  dh  sinh^AdA, 

0  =  Isqs/ns(X0,Yg,zo,A,h,A)[v1a)1+\;2u2]  dh  coshA  sinh  AdA. 
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Generator  Mode  (BNW  Program  Plasma-III).  Critical  scrutiny  of  the 
residual  parametric  integration  structure  of  the  merged  Vlasov-Maxwell 
orbit-integrated  triad  (103-105)  reveals  the  existence  of  one  uniquely 
simple  modal  subgroup,  which  both  admits  full  closed-form  integration  of 
equations  (103-105)  and  yet  contains  just  sufficient  parametric 
flexibility  to  conceivably  be  developable  as  a  generator  function  for  all 
modes.  The  analysis  logic  for  this  generator  mode  (now  operable  on  the 
Cyber-74  computer  as  Battel  1 e-Northwest  Program  Plasma-III)  begins  with 
its  double  delta  function  parametric  ion  distribution, 


ns(WZo’A,’l,,,X)  = 


ns(Xe,Y0,zo,X)  fi(A'-A)  6(h'-h). 


(106) 


For  this  ion  distribution,  the  merged  Vlasov-Maxwell  orbit  integrals 
(103-105)  reduce  to  three  algebraic  equations, 


iK2a2/2u  =  ls  msc2ns(Xe,Ye,z0,X){[P(26)-2h](v1u)1+v2w2) 

4 

-2[v1^(u1)+v2c(m2)]}sinh  A, 

0  =  £sqsns(X0,Ye,zo,A){[d26)-qsb/(Kmsc  sinhA)] (v^+VgWg) 

3 

-w[\)^(w-|)+v2c(u2)]6}sinh  A, 

2 

0  =  Isqsnb(X0,Y0,zo,X)[v1o)1+v2co2]coshA  sinh  A. 


(107) 

(108) 
(109) 


The  first  of  this  algebraic  triad  relates  the  rotary  electromagnetic 
wave  amplitude  to  the  rotary  ion  current  source.  The  last  tv/o  state, 
respectively,  that  both  the  axial  ion  current  and  the  charge  density  are 
zero.  For  a  mixture  of  I  ion  sDecies,  parametric  flexibility  here  spans 
the  real  pair  (a,b)  and  the  41+1  complex  parameters 
[K,(m1,to2,6,A)=(w1,w2,6,A)s],  intricately  interrelated  by  21-1  complex 
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constraints  visible  in  (35-38,50-53)  as  the  I  ion-independent  ratios. 


[qsa/(msc  sinh/Vs)]/[qsb/(Kmsc  sinhAs)]  =  Ka/b, 


(HO) 


and  the  1-1  ion-dependent  ratios, 


Cq-j a/ (m-|C  sinhA]  )]/[q$a/ (m$c  sinhA,.)]  =  (q]ms  sinhA^/fq^  sinhA,), 


(HI) 


to  leave  the  net  equivalent  of  21+3  complex  parametric  degrees  of  freedom 
for  each  2 1 -tuply  indexed  double  period  mode  (v-j.v^)  =  (v-|,v2)s. 

Exploring  this  vast  subset  of  closed-form  transversely-uniform  solutions 
may  take  quite  a  while. 


Three  Ion  Species.  Brief  scrutiny  of  the  three  algebraic  equations 
(107-109)  leads  quickly  to  focusing  initial  investigation  upon  a  mixture 
of  three  ion  species  [such  as  (D+,T+,e  )  for  a  deuterium-tritium  fusion 
plasma,  or  (H  ,D  ,e  )  for  a  planetary  or  stellar  proton  plasma],  and 
simply  inverting  the  3x3  coefficient  matrix  to  obtain  the  parametric  ion 
distribution  set  n  =  (n^n^n^  as 


[mc2{[P(26)-2h](v1aJl+v20J2) 

-1 

iK2a2/2y 

-2[v1?(oj1  )+v2c(w2)]}sinh4A]s 

n(A,h)  = 

[q[[c(26)-qb/(Kmc  sinhA)] (v^+v^) 

• 

0 

> 

-2[v1?(a)^  )+v2;(io2)]<5}sinh3A]s 

k[q(v-|U)^+v2aj2)coshA  sinh^A]^ 

o 

(112) 


where  the  ion-species  index  s  =  (1,2,3)  labels  three  such  column 
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vectors  comprising  the  matrix  to  be  inverted.  The  arguments  (A,Jh)  on  the 
left  stress  the  explicit  appearance  of  the  ion  set  (A,hJ  on  the  right,  as 
distinct  from  implicit  dependence  upon  the  parameter  basis 
(a  ,b ,K,w-j  ,0^2 , 6_)  • 

Transformation  of  this  ion  distribution  from  the  parametric  basis 
back  to  the  wave-frame  (z,A,0,tj;)  basis  involves  now  simply  direct 
substitution  of  the  ion  inner  Hamiltonian  functions  h(z,A,0,ijj)  of  (84) 
for  the  inner  Hamiltonian  parameter  set  h^  on  the  right  side  of  the  last 
equation.  No  density  transformation  determinant  accompanies  such  a  shift 
from  the  parametric  basis  to  the  wave-frame  basis,  since  the  complete 
parametric  solution  (88)  was  defined  in  the  wave  frame.  For  each 
hextuply  indexed  double  period  mode  (vp^),  these  transversely 

uniform  triple  ion  solutions  span  nine  complex  parametric  degrees  of 
freedom,  according  to  the  logic  in  the  vicinity  of  (110,111). 

Four  or  more  ion  species  can  be  handled  quite  similarly,  by  inverting 
any  3x3  coefficient  sub-matrix  appearing  in  the  three  algebraic  equations 
(107-109)  to  express  any  three  ion  distributions  in  terms  of  the  other 
ion  distributions. 


Two  Ion  Species.  Two  ion  species  emerge  from  the  three  algebraic 
equations  (107-109)  as  a  curious  determinantal  degeneracy,  which  can  be 
read  off  the  last  two  homogeneous  equations  (108-109)  by  setting  the 
determinant  of  their  coefficient  matrix  to  zero  to  obtain 


[ { C C ( 26 ) -qb/ ( Kmc  sinhAJKv^+VgU^) 
-2[v-|  c(o>i  )+v2?(aj2)]6}sinhA]s 

[(v^w-i+VgO^JcoshA^ 


013) 
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where  the  ion  species  index  s  =  (1,2)  labels  two  such  column  vectors 
comprising  the  2x2  singular  matrix.  The  parametric  ion  distribution  next 
is  obtainable  by  inverting  the  2x2  coefficient  matrix  for  the  inhomogeneous 
equation  (107)  combined  with  either  of  the  now  coalescent  homogeneous 
equations  (108-109),  to  yield  the  set  n  =  (n^n^  as 

(114) 

where  again  the  ion  species  index  s  =  (1,2)  labels  two  such  column 
vectors  comprising  the  2x2  matrix  to  be  inverted.  For  each  quadruply 
indexed  double  period  mode  (v^  these  transversely  uniform 
double  ion  solutions  span  six  complex  parametric  degrees  of  freedom, 
according  to  the  logic  in  the  vicinity  of  (110,111)  as  modified  to 
include  the  degeneracy  constraint  (113). 
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NUMERICS 

Fully  optimized  digital  computer  generation  of  all  the  preceding 
doubly  periodic  functions  and  their  first  and  second  integrals  proceeds 
[via  Whittaker  and  Watson37^*  by  first  summing  to  solid  (60-bit) 

convergence  thirteen  complex  Fourier  series  defining  the  four  doubly 


quasi-periodic  Jacobi  theta  functions, 

e^z)  =  2£”=q  (-)n  exp[i7Tt(n+l/2)2]  sin[(2n+l )z],  (115) 
62(z)  =  2£”_0  exp[int(n+l/2)2]  cos[(2n+l )z] ,  (116) 
e3(z)  =  1  +  2£”_i  exp[irfn2]  cos[2nz],  (117) 
04(z)  5  1  +  2j“=1  (-)n  exp[iTTfn2]  cos[2ni»],  (118) 

their  four  first  derivatives, 

0j(*)  -  4^=0  (-)n  exp[ iTix (n+1  /2) 2]  (n+1/2)  cos[(2n+l)*],  (119) 
eJU)  =  -4X"=0  exp[inf  (n+l/2)2]  (n+1/2)  sin[  (2n+l  )z],  (120) 
03(z)  =  -4£“_i  expCiirtn2]  n  sin[2nz],  (121) 
84(a)  =  "4In=l  exPti7T,fn2^  n  sin[2nz],  (122) 

and  their  nodal  normal izer  quintet, 

ej(0)  =  4J“=0  (-)"  exp[l7rt(n+l/2)2]  (n+1/2),  (123) 
e'j'(O)  =  -16Co  (_)n  exp[irt(n+l/2)2]  (n+1/2)3,  (124) 
02(0)  -  2^”_q  exp[in1(n+l/2)2],  (125) 
03(0)  =  1  +  2J“„1  exp[iirtn2],  (126) 
0,(0)  =  1  +  2£„,  (-)"  exptliftn2].  (127) 


Here  the  primitive  (tiniest)  complex  periods  (2w-|,2u)2)  are  to  be  ordered 
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in  phase  angle  (by  relabeling,  or  by  replacing  either  period  with  its 
negative)  so  that  the  primitive  period  ratio  f  i  o>2/u>.|  has  a  positive 
imaginary  part  to  provide  uniform  series  convergence.  Because  of  the 

o 

presence  of  the  double-exponential  attenuation  factor  exp[iirfn  ],  all 
these  series  display  ultra-rapid  convergence  over  the  full  complex  field, 
except  quite  near  their  dense  line  of  intrinsic  essential  singularities 
along  Im(t)  -  0. 

The  doubly  periodic  Weierstrass  pay  function  P(z),  the  additive 
doubly  quasi-periodic  Weierstrass  zeta  function  ;(i),  and  the 
multiplicative  doubly  quasi-periodic  Weierstrass  sigma  function  o(z)  all 
then  are  generated  quite  simply  by  a  flexible  variety  of  sequences 


typified  by  [but  by  no  means  limited  to39^Sec* 

e}  =  P(a)1)  =  (tt/2^  )2[03(O)+eJ(O)]/3,  (128) 

e2  =  P(w2)  =  -(^^^(oH^O)]^,  (129) 

e3  =  P(u3)  =  (TT/2u)1)2[e2(0)-eJ(0)]/3,  (130) 

n,  =  d^)  =  -(TT/2)2[e“  ’  (0)/3e  ■  (o)3/u>1 ,  (131) 

n2  =  dw2)  =  n-|^-i7r/2o>1 ,  (132) 

n3  =  c(w3)  =  -n1-n2»  (133) 

P(h)  =  e1+[(Tr/2u)1)e2(TTH/2u).|  )ej(0)/e1(Tra/2a)1)e2(0)]2,  (134) 

P(z)  =  e2+[(TT/2aJl)e4(TTZ/2a)1)ej(O)/01(7ra/2a)1)e4(O)]2,  (135) 

P(z)  =  e3+[(7T/2u)1)e3(TTZ/2u,1)0](O)/e1(TTZ/2a)1)03(O)]2,  (136) 

c(^)  =  ni2/u)i  +  (TT/2w1  )0]  (ttz/2o)1  )/©i  (Trz/2a)i ) ,  (137) 

a(z)  =  (2u31/tt)[01(ttz/2uJi)/0](O)]  exp[n122/2a>1].  (138) 
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Logarithmic  differentiation  of  the  transparent  identity, 

p'(z)  =  -2(71/2^  )3[e](O)]2[02(TTZ/2u)1)03(TTZ/2w1)e4(TTZ/2u1)/0^(TT2/2u1)],  (139) 

followed  by  cyclic  half-period  phase  shifts,  generates  the  quartic  roots 
(50-53)  from  the  symmetric  sequence, 

P"  (z)/P '  (2)  =  (tt/2u)-|  )[02(tt2/2cj1  )/02(tiz/2o)1  )+0^(ttz/2u1  )/03(ttz/2u1  ) 

+e^(irz/2o)1  )/04(ttz/2u)1  )-30j  (ttz/2w-|  )/0^  (ttz/2u1  )] ,  (140) 

P"  (z+u1  )/P'  (z+w1 )  =  (tt/2(a)-j  ) [©3 (ttz/2oj-j  )/03(ttz/2u»1  )+04(7iz/2oj1  )/04(ik  / 2 ) 

+G]  (nz/2u)1  )/01  (ttz/2oi1  )-30J(ttz/2lo1  )/02(uz/2(Jj1  )] ,  (141 ) 

P"(z+u)2)/P'  (z+co2)  =  (7t/2o>1  )[6| (ttz/2co1  )/01  (ttz/2u)1  )+0^(tth/2u)-1  )/02(ttz/2uj1  ) 

+0^(n2/2a)1  )/03(ttz/2co1  ) -30^712/2^  )/04 (ttZ/203-,  )] ,  (142) 

P"  (z+u3)/P'  (z+ui3)  =  (tt/2o)-|  )[G4(7t2/2wj  )/04(tt2/2w1  )+0]  (t7z/2u)1  )/@1  (t7Z/2o)1  ) 

+0  2  (  ttz/  2lo  1  )/02(7T2/2u)1  )-30^(ttz/2u)1  )/03(t72/2u1  )] .  (143) 

All  these  numerics,  it  must  be  stressed,  are  based  directly  upon  the 
cyclicly  symmetric  Weierstrass  definitions  (as  developed  most  recently  by 
Whittaker  ),  both  to  minimize  analysis  labor  and  to  simplify  digital 
computer  loops.  [The  cyclicly  asymmetric  Jacobi  definitions  (as 
developed  most  recently  by  Southard39)  cannot  survive  optimality  scrutiny 
with  respect  to  either  analysis  labor  or  digital  computer  efficiency.] 


r 
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PARAMETER  VARIATION 

Variation  of  parameters  [such  as  the  parametric  basis 
(A,a ,b,K,aji  .ajg ,6.)  spanning  the  preceding  Vlasov-Maxwell  ion-plasma 
solutions]  to  meet  intricately  coupled  constraints  [such  as  the  cascaded 
horrors  (35-38,50-53,110-111,113)]  involves,  in  general,  varying  a  subset  a 

of  input  parameters  to  obtain  a  subset  J3(a)  of  output  responses  which 
satisfy  a  subset  of  constraints  x(B)  =  0.  Mathematically,  this  means 
inverting  the  cascaded  vector  function  equation, 

l(£(a))  =  0,  (144) 

to  obtain  a.  In  the  context  of  machine  analysis,  a  consists  of  some  of 
the  input  data,  8(a)  consists  of  some  of  the  output,  and  x(£)  is  a 
specified  evaluation  function  measuring  the  extent  to  which  the 
parameters  a  depart  from  constraint  requirements. 

Consider  looping  the  entire  machine  analysis,  by  varying  a 
systematically  in  successive  passes  by  small  input-controlled  increments 
A(a)  to  generate  a  matrix  difference-ratio  representation  of  the 
error/value  sensitivity  gradient  [A(a)]”1 '[A(x)] »  and  then  expanding  the 
constraints  (144)  in  the  form 

l(B(u))  +  (da)*[A(a)]_1 *[A(x)]  =  0.  (145 ) 

Solving  this  set  of  locally  linear  algebraic  equations,  by  cascaded 
matrix  inversion,  yields  for  the  parameter  changes  da  needed  to  meet 
constraints, 

da  -  -  x(B(a))-[A(x)]-1 *[A(a)].  (146) 

It  is  to  be  noted  here  that  the  parameter  vector  a  and  the  constraint 
vector  x  must  have  the  same  dimensionality  and  that  the  intermediate 
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vector  6  must  have  dimensionality  equal  to  or  greater  than  that  common 
to  a  and  x*  Applying  the  parameter  refinement  da  to  the  parameter  basis 
a  now  generates  the  recursive  refinement  algorithm  for  the  variational 
logic, 


a  =  a  - 


(147) 
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CONCLUSIONS 

Exact  expressions  for  nonlinear  ion  orbits  under  the  influence  of  an 
axial  magnetic  field  and  an  arbitrary  rotary  wave  have  been  derived  in 
terms  of  complex  elliptic  functions.  As  a  direct  byproduct  of  these 
momentum  and  position  space  expressions,  the  complete  set  of  six  non¬ 
trivial  constants  of  the  motion  were  identified.  The  constants  permitted 
the  Vlasov  equation  to  be  satisfied  identically  by  a  flexible  phase-space 
distribution  function  of  the  constants.  This  functional  form  allowed 

Maxwell 's  equations  to  be  integrated  completely  to  obtain  a  nontrivial 
eigenvalue  equation. 

Exact  nonlinear  expressions  for  the  phase-space  distribution  are 
obtained  explicitly  for  two-  and  three-ion  species  in  terms  of  six  complex 
eigenvalues.  Variation  of  startinq  values  that  will  yield  physically 
acceptable  eigenvalues  is  not  trivial  since  their  determination  must  be 
evaluated  over  a  set  of  multiply  indexed  double  period  modes.  The  evalua¬ 
tion  of  the  complex  elliptic  functions  is  facilitated  by  the  series 
expansions  which  converge  rapidly.  The  numerical  evaluation  of  the  proper 
Physical  modes  remains  a  difficult  problem,  but  an  exhaustive  computer 
program  has  been  formulated.  The  nature  of  the  problem  necessitates  a 
parameter  variation  procedure  in  the  complex  eigenvalue  hyperspace. 

Several  applications  of  these  results  remain  to  be  explored.  Various 
techniques  for  stimulation  of  UIF-VLF  (0.1-105  Hz)  radio  noise  in  the 
magnetosphere  and  ionosphere  require  a  nonlinear  description.  Hhen  specific 
experiments  are  contemplated,  a  thorough  examination  of  the  adjacent 
parameter  domain  would  be  advisable  using  the  foregoing  nonlinear  solutions 

in  order  to  determine  the  optimum  natural  conditions  for  good  diagnostic 
signatures. 


» 
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Fusion  research  using  dense  plasmas  may  also  be  served  in  various 
ways.  The  existing  formalism  might  be  useful  for  optimization  of  condi¬ 
tions  for  electromagnetic  pumping  of  plasma  energy  density.  Radiation 
loss  rates  for  simple  geometries  can  be  evaluated  at  their  nonlinear 
limit.  Finally,  the  nonlinear  theory  may  provide  improved  interpretation 
of  observed  radio  noise  emissions  which  could  then  be  developed  as  a 
diagnostic  tool. 
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